Heat capacity of matter beyond the Dulong-Petit value 
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We propose a new simple way to evaluate the effect of anharmonicity on a system's thermodynamic 
functions such as heat capacity. In this approach, the contribution of all potentially complicated 
anharmonic effects to constant-volume heat capacity is evaluated by one parameter only, the coeffi- 
cient of thermal expansion. Importantly, this approach is applicable not only to crystals but also to 
glasses and viscous liquids. To support this proposal, we perform molecular dynamics simulations 
of several crystalline and amorphous solids as well as liquids, and find a good agreement between re- 
sults from theory and simulations. We observe an interesting non-monotonic behavior of liquid heat 
capacity with a maximum, and explain this effect as a result of competition between anharmonicity 
at low temperature and decreasing number of transverse modes at high temperature. 



INTRODUCTION 

One of the central and most recognizable results of 
statistical physics is the value of constant-volume heat 
capacity, C„, of a harmonic and classical solid: 



(1) 



where N is the number of atoms and fee = 1. Known 
as the Dulong-Petit law, Eq. (1) is the result of a solid 
having phonons I]. 

Experimentally, Cy is almost never 3A^ even in the clas- 
sical limit 1, where wd is Debye frequency, an 
effect attributed to anharmonicity of interatomic interac- 
tions In addition to heat capacity, anharmonicity 
governs many other properties of condensed matter sys- 
tems, including thermal expansion, thermal and electric 
conductivity, elasticity, phase transitions, defect mobil- 
ity, melting and so on. 

There has been a large amount of research into anhar- 
monic effects 0, 3, 0| that has resulted in qualitative un- 
derstanding of the effect anharmonicity on system prop- 
erties. The common approach is to expand the potential 
energy U in Taylor series over atomic displacements u: 
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where the anharmonic coefficients (fix... are given by the 
derivatives at equilibrium separations in a usual way fsl] . 



As noted by Cowley 0, 4>x... are very complicated to 
evaluate even if the potential functions are known. Com- 
plications related to evaluating 4>x... necessitated approx- 
imations, which, as Cowley further notes 0, are quite 
inadequate for real systems and are useful in order-of- 
magnitude calculations only. However, assuming that 
interactions include only pair and short-range (nearest- 
neighbor) interactions (/> and considering, for example, a 
face-centered cubic lattice, low-order perturbation theory 
gives Cy as a function of (f) and T as 
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This relationship is one of the few that provide a closed 
form for evaluation of C^, assuming that (fi are known 
and, importantly, represent a faithful representation of 
interatomic forces. 

Unfortunately, the quantitative evaluation of anhar- 
monicity effects has remained a challenge, with the fre- 
quent result that the accuracy of leading-order anhar- 
monic perturbation theory is unknown and the mag- 
nitude of anharmonic terms is challenging to justify 
[3, H, H, [nl Experimental data such as phonon 

lifetimes and frequency shifts can provide quantitative 
estimates for anharmonicity effects and anharmonic ex- 
pansion coefficients in particular, although this involves 
complications, and limits the predictive power of the the- 
ory Q . As noted starting from the early studies [1, [gl , 
the main problem with the approach based on expansions 
such as Eq. (2) and subsequent understanding of anhar- 
monic effects is that good-quality models for interatomic 
forces are not generally available. 

The problem of the anharmonic theory relying on the 
knowledge of interatomic interaction models has been 
noted earlier 0, [ll[ . It has been stated that "undoubt- 
edly the unsatisfactory nature of these models is the lim- 
iting factor in our understanding of many anharmonic 
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properties", and that if anharmonic calculations are to 
have any quantitative significance, the realistic models 
are necessary Consequently, theoretical work on in- 
teratomic potentials was stated as an essential future ef- 
fort at the time We note in passing that despite the 
progress in materials modeling since that time, the prob- 
lem remains. Indeed, apart from relatively small number 
of materials, a negligible fraction of all known ones, it 
has proved impossible to develop a general recipe for suc- 
cessfully mapping the interatomic interactions onto the 
sets of empirical functions to be used in expansion such 
as (2). The problem is particularly acute with modern 
materials which often have complicated interactions in 
the form of hydrogen-bonded, polymeric and many-body 
interactions, magnetic correlations, non-trivial band gap 
changes with temperature, anisotropy, layered structures, 
large number of distinct atoms in organic and biological 
systems and other factors that severely limit the devel- 
opment of high-quality interaction models. Importantly, 
because the anharmonic effects are believed to be small, a 
small departure of a potential model from a high quality 
one renders the partitioning into harmonic and anhar- 
monic parts (2) and subsequent interpretation of anhar- 
monicity meaningless. 

Another approach to treat anharmonicity is to in- 
voke Griineisen approximation, where the softening of 
phonon frequencies uJi is quantified by parameters — 
~'uj (w')t' ^^'^ discuss the macroscopic equations of 
state However, Cy was not previously calculated in 
this approach in the form free of adjustable parameters 
and suitable for direct numerical evaluations. 

In view of persisting difficulties of evaluating anhar- 
monic effects, it is important to have an alternative gen- 
eral method of estimation of anharmonic Cy. It is also 
important to have a method applicable not only to crys- 
talline systems, but also to amorphous solids as well as 
liquids, systems for which the traditional perturbation 
approaches are not suitable, as discussed below in more 
detail. 

In the course of studying the problem of glass tran- 
sition, we have proposed [IJ] that the effects of anhar- 
monicity on Cy can be evaluated as 

Cy = 3N{1 + aT) (4) 

where a is the coefficient of thermal expansion. 

There is no contradiction in the relationship (j4|) be- 
tween the constant-volume Cy and thermal expansion, 
as might be perceived. As discussed below in detail, the 
relationship is due to the softening of bulk modulus with 
temperature at constant volume due to intrinsic anhar- 
monicity, an effect that can be related to a in Griineisen 
approximation. 

In Eq. ([4]) , all potentially complicated effects of anhar- 
monicity discussed above are evaluated by one parame- 
ter, a. Importantly, a is not an adjustable parameter. 



but is fixed by system properties. Another important 
feature of Eq. (|4]) is that a can be independently mea- 
sured or calculated in a straightforward way no matter 
how complicated interactions in a system are. Appeal- 
ingly simple, Eq. ^ provides an important and straight- 
forward way of estimating the effect of anharmonicity on 
Cy . Perhaps not unexpectedly, the simplicity is achieved 
by making approximations, and this paper is partly de- 
voted to assessing these approximations, a point to which 
we return below. 

Importantly, Eq. (|H) can be used to evaluate anhar- 
monicity in two important types of condensed matter sys- 
tems, glasses and liquids, for which calculations based on 
anharmonic expansions such as (2) do not work. Indeed, 
the evaluation of the anharmonic terms in Eq. (2) and 
coefficients ^x... involves sums over wave vectors fc in a 
crystal 0, y] . On the other hand, k are not defined in 
amorphous glasses, at least not at large k. In liquids, an 
expansion such as (2) can not be made even in principle 
because atoms do not oscillate around fixed positions as 
in solids, which is the starting point of theories based on 
Eq. (2) and similar ones. 

In this paper, we extend our new approach, and ad- 
dress the validity of Eq. (|4]) across a wide range of crys- 
tals, glasses and viscous liquids. We perform molecular 
dynamics (MD) simulations, and find a good agreement 
between simulation results and Eq. @ in several crys- 
talline and amorphous solids as well viscous liquids in a 
wide temperature range. 

We note that using MD simulation to study Eq. (|4]) 
has two important advantages over experiments. First, 
experimental Cy is calculated from the measured Cp as 
Cy = Cp~ VTa^B, where B is the bulk modulus. There 
are uncertainties in experimentally determined a and B, 
particularly at high temperature, which implies uncer- 
tainty in C„ [1]. In the MD simulation, this problem 
does not originate because simulations can be performed 
at constant volume. Second, the classical limit <^ 1 
giving Cy — SN is not achieved in many experimental 
systems due to high cjd [3] . Consequently, it is often not 
clear to what extent the deviation of experimental Cy 
from 3N is due to anharmonicity or quantum effect of 
phonon excitation. This issue does not originate in our 
MD simulations, which are classical. 

We finally note that when evaluations of anharmonic 
effects are possible for certain systems, traditional per- 
turbation approaches achieve the accuracy at the level of 
order-of-magnitude agreement with experiments or sim- 
ulations (see, e.g., Refs 0, H, [Hj). We aim for at least 
the same level of accuracy in our new general method 
of evaluating anharmonic effects. The accuracy is deter- 
mined by certain approximations that are used to derive 
a simple form of Eq. 1^. We find that Eq. gives 
correct order-of-magnitude evaluation of anharmonic ef- 
fects, the result that is considered as best of what can be 
achieved in the traditional perturbation expansion ap- 
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proximations. 



THEORY 



We start with the derivation of Eq. The free 

energy of a harmonic soUd in the high-temperature ap- 
proximation is F — SA^Tln where u)^^ — a^W2---W3Af 
is geometrically averaged phonon frequency ,1.]. In the 
harmonic case, w is constant, giving the entropy S = 



(9F) 



■iNn+ln-^) and 



Anharmonicity results in the decrease of oj with temper- 
ature. Then, S* = 37V (l-f In " 



i3f).and 
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where the derivatives are taken at constant volume. 

In the high-temperature limit where F = 3NT In ^ , 
Eq. ([5]) is exact, and is the starting point of our theory. 
Evaluation of Cy requires the knowledge of ^ , which we 
calculate below. 

The phonon pressure, Pph, is Pph - (|f = 2^^, 

3N 

where 7 is the average Griineisen parameter 7 = ^ 7^ 

1=1 

and "fi = —-^ (w")t ^^^^ gives the bulk modulus 
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where 



V "^^^ \^ dT 

q = gi^v • Experimentally, g is known to be fairly con- 



stant across the range of systems (e.g. q=2.l for Pb, 3.2 
for Ge 0, 1.4 for MgO [J, 1.5-2 for alkali hahdes 
1.7 for MgSiOs perovskite fl6^ etc). For simplicity, we 

set (^^§jr-^ — as this does not affect our order-of- 

magnitude evaluations of c^, a point to which we return 
below. Using 7 = ^^^^ and B = Bq + i?ph, where B and 
Bq is the total and static bulk modulus, respectively, we 
find 



dB 



ph 



dT 



-a{Bo + Bph) 



(6) 



where we set Cy = 3A^ in this approximation. 

For small aT, which is often the case in the exper- 
imental temperature range, Eq. JE]) implies B cx — T, 
consistent with the experiments [4|]. We note that ex- 
perimentally, B linearly decreases with temperature at 
both constant pressure and constant volume (constant- 
volume decrease can be small in some systems) [i,[i3-[il- 
The decrease of B with T at constant volume is due to 
the intrinsic anharmonicity related to the softening of in- 
teratomic potential at large vibrational amplitudes; the 
decrease of B at constant pressure has an additional con- 
tribution from thermal expansion. 

The next step is to assume that cx; B, a relationship 



because ojf = k'^(? oc B -I- |G and the shear modulus G 
scales with B via the Poisson ratio that is nearly constant 
in all systems. Therefore, cx B is applicable to any 
system as long as the phonon spectrum is treated in De- 
bye approximation, as is often the case. In a general case 
of a spectrum that includes optic modes, the relation- 
ship up oi B can be addressed by studying how uji and B 
change in response of external parameters such as tem- 
perature and pressure. It has been found that cuf (x B 
is the case for optic modes in a wide temperature range, 
both longitudinal and transverse The increase of 

uji including acoustic and optic modes is also seen in a 
wide pressure range, accompanied by the simultaneous 
increase of B |2l| . 

Finally, combining uj'^ Bq + -Bpii and Eq. (j6|), we find 
i (^)^ = -|. Putting the last relationship in Eq. © 
gives Eq. (|4]). We note that the last two terms in Eq. 
© cancel out if (jf )^ oc w, as is the case here. 

As follows from the previous discussion, the evaluation 
of Cy can be made more precise if values of q are retained 

in the calculation. In this case, (^^§p-^ = — (5(i3o+Sph), 

where 6 = a{q — 1). Gombining it with oj'^ (x Bq + _Bph 
gives 4 (^)^ = — |. Using it in Eq. ([5]) gives 



Cy = 3A^(1 + ST) 



(7) 



Here, similar to Eq. ([4]), all anharmonic effects are rep- 
resented by one parameter, 6. This parameter quantifies 
the decrease of B with temperature at constant volume. 
Goncerned with demonstrating an order-of-magnitude 
evaluation of Cy using our new approach, we will not pur- 
sue Eq. dl]) further, and concentrate on Eq. (jl]). 



MOLECULAR DYNAMICS SIMULATIONS 

We now discuss our MD simulations. We have aimed 
for diversity of structures and interactions, and conse- 
quently chosen several systems with different symme- 
try, structure and interatomic potentials: crystalline Ge, 
NaGl, AI2O3 (corundum), TiOs (rutile), ZrSi04, SiOs 
glass and a model liquid system. For Ge, we used many- 
body environment-dependent ( "bond-order" ) Tersoff po- 
tential [2^. Here, the interaction strength between any 
two atoms depends on their environment and coordina- 
tion. This potential therefore represents an example of a 
crystalline system which can not be treated in the expan- 
sion approach (2). Empirical potentials for AI2O3 124 



23-ii| and Si02 glass [3C 



Ti02 |25l, NaGl [261] ZrSi04 
included long-range Goloumb and short-range Bucking- 
ham or Morse interactions. Ref. jsij discusses details of 
generation of Si02 glass structure. For the liquid, we em- 
ployed Lennard- Jones (LJ) potentials designed to simu- 
late a binary liquid in the supercooled viscous state [3^ . 



that holds true if ujf oc B. For acoustic modes, uf oc B The binary liquid consists of two distinct atomic types 
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with different interaction parameters and effective sizes 
to avoid crystallization at low temperature. 

We note here that the empirical potentials we em- 
ployed may or may not closely reproduce the experimen- 
tal a or other properties such as c-a or B. However, this 
is not important for our study as we aim to show that a 
given force field, even though approximate, still results 
in the relationship between c„ and a given by Eq. (|4]). 
In this sense, it is only important that a force field gives 
physically sensible set of Wi (e.g., real w^) and other phys- 
ical characteristics such as elasticity and thermal expan- 
sion that show commonly observed temperature depen- 
dence, because our derivation of Eq. ([4]) is based on these 
properties and relationships between them. 

We have used DL_POLY programme [s^ for our MD 
simulations. For solids, the number of atoms was in 
the 12,000-27,000 range depending on the system. For 
the LJ liquid, we used 64,000 atoms. We have verified 
that increasing the number of atoms does not change the 
results. The energy of the system, E, was calculated 
in constant-energy and volume ensemble simulations by 
equilibrating the system at a given temperature. The sys- 
tem volume and a were calculated in constant-pressure 
ensemble simulations. We have performed simulations in 
a wide temperature range (see Figure 1) with tempera- 
ture step of 1 K. Each temperature point was simulated 
on a separate processor using our high-throughput com- 
puting cluster. The constant-volume specific heat was 
calculated as c„ = To reduce the fluctuations of 

the derivative, we have fitted the energy using high-order 
polynomials and cubic splines, and verified that c^, is not 
sensitive to the polynomial order used and fitting param- 
eters. 

In Figures 1-2 we show the calculated c„ and relative 
volume where Vq is the system volume at the lowest 
simulated temperature, for 6 solid systems and for the LJ 
liquid. We observe that for different systems increase 
above the Dulong-Petit value of 3 in the wide tempera- 
ture range. We found an exception to this behavior in 
crystalline Ar, where the increase of c« is preceded by 
its decrease at low temperature. In soft crystals such as 
Ar with large anharmonicity (7 =3.5), we do not expect 
the starting assumptions of the Griineisen approximation 
that we employed here to hold. Figures 1-2 enable us to 
see how well the slope of predicted by Eq. (jj]) agrees 
with the actual value of a. Consequently, we calculated 
etc from Figure la as = 3(1 -I- ttcT) and a = -^-^y 
from Figure lb and 2b. For the LJ liquid, etc was calcu- 
lated from the linear increase of at low temperature 
in Fig. 2a, for the reasons discussed below in detail. For 
some systems, and are not linear with tempera- 
ture in the whole temperature range. In this case, we 
have calculated ac and a at each temperature, and have 
taken the average. 

The calculated values of Uc and a are: crystalline Ge 
(ac = 3.6 • 10-5 K-i, a = 2.6 • 10"^ K"!), TiOa (ac = 
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FIG. 1: Cv (a) and (b) for simulated crystalline and amor- 
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1.1 • 10-5 K-i, a = 2.8 • 10-5 K-i), NaCI (o^c = 7 • 10" 



14 • 10-5 K-i), ZrSi04 (ac 
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1 --0.7- 

5 



a : 



a = 2 • 10-5 K-i), AI2O3 (ac = 1-3 • 10-5 K 
10-5 K-i), SiOz glass (ac = 2.4-10-5 K-^, a = 2.9-10 
K-i), LJ liquid (ac = 1.75 - IQ-''' R-^. a = 1.72 - IQ-^ 
K-i). 



We observe that the anharmonic contribution to can 



be evaluated by Eq. ([4]) fairly well. 



, averaged over 



all systems, is within 40%. The discrepancy between the 
predicted and calculated values is within the approxima- 
tions we made in deriving Eq. ^ , including our neglect- 
ing the dependence of 7 on volume, which can alter the 
predicted ac by up to a factor of 2 (see Eq. ([7]) and the 
discussion preceding Eq. ([6|)). We therefore find that Eq. 
(HI) gives the correct order-of-magnitude evaluation of an- 
harmonic effects, the result that is considered as best of 
what can be achieved using the traditionalperturbation 
expansion approximations (see, e.g., Refs [1^0, [HI). 
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FIG. 2: Ci, (a) and ^ (b) for LJ liquid, (c) shows coordinates 
of three atoms with large atomic displacements. 



HEAT CAPACITY OF A LIQUID IN THE 
VISCOUS REGIME 



We now discuss why and how the above theory ap- 
plies to viscous liquids in addition to solids. We first 
demonstrate that the temperature range where lin- 
early increases in Fig. 2a, and where we calculate ac, 
corresponds to a viscous liquid. Following a somewhat 
general definition, a viscous liquid is a liquid whose re- 
laxation time r is much larger than Debye vibration pe- 
riod of about 0.1 ps: r ;» td. Here, r is the average 
time between consecutive atomic jumps in a liquid at one 
point in space In Fig. 2c we plot the coordinates of 
three atoms from the simulation of the binary LJ liquid 
at 50 K, corresponding to temperature in the middle of 
the linear increase of Cy in Fig. 2a. We observe atomic 
displacements reaching 8 A during the time of our simu- 
lation, witnessing that large- amplitude diffusive motions 
are present as is the case for liquids. Second, we estimate 
T as the average time between atomic jumps by defini- 
tion. Averaged over different atoms, t is approximately 
15 ps. Therefore, r td, corresponding to a viscous 



liquid. 

There are two types of motion in a liquid: phonon mo- 
tion that includes one longitudinal mode and two trans- 
verse modes with frequency cj > i, and diffusional mo- 



tion [22]. Consequently, the total liquid energy, E, is 
the sum of the phonon energy, -Bph, and diffusional en- 
ergy, -Edif : E = i?ph -I- Edii- Edif includes both kinetic 
energy of diffusing atoms and potential energy of their 
interaction with other atoms. As argued by Frenkel, a 
particle spends time r vibrating in between jumps 22 1. 
The time it takes a particle to jump from one equilibrium 
position to the next is approximately equal to td . There- 
fore, the probability of a jump is p — In statistical 
equilibrium, the number of atoms in the transitory dif- 
fusing state is A^'dif — Np, where N is the total number 
of atoms, giving 



iVdit = N— 

T 



(8) 



Eq. ([5]) implies that in a viscous liquid where r ^ 
Td, the relative number of diffusing atoms at any given 
moment of time is negligible. Consequently, -Edif can be 
ignored, giving E = i?ph at any given moment of time. 
It is easy to show that the same result, E = Ep^, also 
applies to the energy averaged over time r 

The phonon energy of a liquid in the regime r 3> 
is given, to a very good approximation, by the phonon 
energy of its solid. This is supported by the explicit equa- 
tion for the liquid energy in the next section, and can be 
qualitatively discussed as follows. The only difference be- 
tween the phonon states in a liquid and a solid is that 
the former does not support all transverse modes as a 
solid does, but only modes with frequency u! > [2^ . 
When T Td , the fraction of missing transverse modes 
in a liquid is negligible and, furthermore, contributes a 
vanishingly small term to the phonon energy because the 
phonon density of states is proportional to w^. 

We therefore conclude that the energy of the viscous 
liquid is equal to the phonon energy, as in the solid. 
Consequently, Eq. derived for solids on the basis 

of phonons and Griineisen approximation, applies to vis- 
cous liquids too. This explains our earlier finding that 
the increase of liquid Cy in the low-temperature viscous 
regime in Fig. 2a is well described by our proposed Eq. 
©• 



THE ORIGIN OF NON-MONOTONIC 
BEHAVIOR OF LIQUID c„ 

It is interesting to note the non-monotonic behavior of 
Cv in Fig. 2 with a maximum. We explain this behav- 
ior as a result of two competing effects. On one hand, 
Cy increases in the viscous regime due to anharmonicity 
as discussed above. On the other hand, Cy decreases at 
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high temperature as a result of progressively decreasing 
number of transverse waves with frequency oj > -. We 



have studied this effect in a series of recent papers [35- 
37[, and shown that the associated decrease of Cy is in 



quantitative agreement with experimental data of many 
liquids. Explicitly, the energy of a classical liquid is |36j : 



be used to evaluate anharmonic in a system of any 
complexity including glasses and viscous liquids, in con- 
trast to previous treatments of anharmonicity. In liquids, 
anharmonicity results in the increase of Cy at low tem- 
perature, contributing to a non-monotonic behavior and 
a maximum of c„. 



E = NT 1 



aT 



(9) 



At low temperature when r ;» td, Eq. © gives 
E = 3iVr(l ^) and = 3iV(l + aT), Eq. (4). 
This is the result we observe in Fig. 2a at low tempera- 
ture. At high temperature when r — s> td, the last term 
in Eq. (|9]), ^3 — (^)'^^ , can not be ignored. Its decrease 

with temperature dominates over T (l -I- ^) because r 
decreases with temperature exponentially or faster. The 
result is that in the low- viscous regime r — td, de- 
creases with temperature [35l - l37l |. The combination of 
two competing effects gives the maximum of c„ as is seen 
in Fig. 2a. 

We finally note that experimentally, the non- 
monotonic behavior of Cv shown in Figure 2a is challeng- 
ing to observe. On one hand, the noticeable decrease of c« 
requires r approaching td as discussed above and, there- 
fore, requires experimenting with low- viscous liquids such 
as metallic, noble-atom and some molecular liquids 37 1. 
These liquids tend to easily crystallize on cooling, pre- 
venting the formation of the viscous regime and accom- 
panied linear increase of c„. On the other hand, the 
linear increase of c„ could be observed in viscous liquids 
such as silicates, chalcogenide and other systems. How- 
ever, reaching low- viscous regime r — td in these sys- 
tems and accompanied decrease of Ct, due to the loss of 
transverse waves requires high temperatures where ex- 
periments are challenging. Moreover, viscous liquids of- 
ten have strong bonds and high Debye temperature of 
internal vibrations, with the result that continues to 
increase even at high temperature due to progressive ex- 
citation of internal vibrations, counteracting the decrease 
of Cy due to the loss of transverse modes. As a result, ex- 
periments typically observe either decrease of Cy in low- 
viscous liquids or increase of c„ in high-viscous liquids 
but not both. These problems did not originate in our 
MD simulations in which we were able to reach both low- 
viscous (r — >• Td) and high- viscous liquid state (r td) 
and which, furthermore, were classical. 



SUMMARY 

In summary, we have discussed a new way of evaluating 
the effects of anharmonicity on system's thermodynamic 
functions such as heat capacity, and have demonstrated 
its good predictive power. Importantly, our theory can 
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